home *** CD-ROM | disk | FTP | other *** search
/ Our Solar System / Our Solar System.iso / miscprog / astrom / planets.pas < prev    next >
Pascal/Delphi Source File  |  1985-02-23  |  29KB  |  1,001 lines
  1. PROGRAM PLANETS;
  2. {
  3.  
  4.                                Version 1.0
  5.  
  6.      This Program computes information relating to the position, distance,
  7.   magnitude, etc for the major planets on a specified date and time.
  8.   Refer to PRACTICAL ASTRONOMY WITH YOUR CALCULATOR by Peter Duffett-Smith
  9.   for the calculation methods.
  10.  
  11.      This program is placed in the public domain "as is".
  12.  
  13.      Comments may be sent to the author:
  14.  
  15.              Larry Puhl
  16.              6 Plum Court
  17.              Sleepy Hollow, Ill. 60118
  18. }
  19.  
  20. CONST
  21.    W = 12; {Width of real printouts}
  22.    D = 6;  {Decinal places in real printouts}
  23.    Number_of_Planets = 9;
  24.    CenterX = 300;
  25.    CenterY = 92;
  26.    Screen_Aspect = 2.7;
  27.    Scale_for_inner_planets = 60;
  28.    Scale_for_outer_planets = 2;
  29.  
  30.    Display_Degrees_in_decimal : boolean = false;
  31.    Daylight_Savings_Correction_Enabled : boolean = true;
  32.    Display_Hour_Angle_in_Degrees : boolean = false;
  33.    Display_RA_in_Degrees : boolean = false;
  34.    Zone_correction : integer = 6;
  35.    Longitude : real = 88.3;
  36.    Latitude : real = 43.1;
  37.  
  38.  
  39. TYPE
  40.    orbit =
  41.       RECORD
  42.          name : STRING[10];
  43.          Epoch : integer;
  44.          Period : real;
  45.          Long_at_Epoch : real;
  46.          Long_of_Per : real;
  47.          Ecc : real;
  48.          Axis : real;
  49.          Inc : real;
  50.          Long_of_Ascend_Node : real;
  51.          Size : real;
  52.          Brightness : real;
  53.       END;
  54.    MSDOS_REGISTERS =
  55.       RECORD
  56.          AX,BX,CX,DX,BP,SI,DI,DS,ES,FLAGS: Integer;
  57.       END;
  58.  
  59. VAR
  60.    Bios_Display_Mode : Integer Absolute $0000 : $0465;
  61.    MSDOS_time : MSDOS_REGISTERS;
  62.    ch : char;
  63.    zero : real; {used to print 0}
  64.    planets : ARRAY [1 .. 10] OF orbit;
  65.    num : integer;
  66.    day,month,year,day_of_year : integer;
  67.    JD : real;
  68.  
  69.  
  70.    planet_long_of_ascend_node : real;
  71.    planet_mean_anomaly : real;
  72.    planet_heliocentric_long : real;
  73.    planet_true_anomaly : real;
  74.    planet_radius_vector : real;
  75.    planet_helio_ecliptic_lat : real;
  76.  
  77.    earth_long_of_ascend_node : real;
  78.    earth_mean_anomaly : real;
  79.    earth_heliocentric_long : real;
  80.    earth_true_anomaly : real;
  81.    earth_radius_vector : real;
  82.  
  83.    GMT : real; {Greenwich mean time}
  84.    GST : real; {Greenwich Siderial time}
  85.    LCT : real; {Local Civil time}
  86.    LST : real; {Local Siderial time}
  87.  
  88.    Hour : integer;
  89.    Minute : integer;
  90.    T,B : real;
  91.  
  92.    projected_heliocentric_long : real;
  93.    projected_radius_vector : real;
  94.  
  95.    Geocentric_ecliptic_Long : real;
  96.    Geocentric_latitude : real;
  97.  
  98.    RA : real;{equatorial coordinates}
  99.    DEC : real;
  100.  
  101.    Phase : real;
  102.    Distance_from_Earth : real;
  103.    Diameter : real;
  104.    Magnitude : real;
  105.  
  106.    Hour_Angle : real;
  107.    Hour_Angle_in_Degrees : real;
  108.    Azimuth : real;
  109.    Altitude : real;
  110.  
  111.    LST_Rise : real;
  112.    LST_Set  : real;
  113.    GST_Rise : real;
  114.    GST_Set  : real;
  115.    GMT_Rise : real;
  116.    GMT_Set  : real;
  117.    LCT_Rise : real;
  118.    LCT_Set  : real;
  119.    Azimuth_Set : real;
  120.    Azimuth_Rise : real;
  121.  
  122.    y : real; {intermediate variable}
  123.    x : real; {intermediate variable}
  124.    z : real; {intermediate variable}
  125.    zmax : real; {intermediate variable}
  126.    time : real;
  127.    Daylight_savings_start : real;
  128.    Daylight_savings_end  : real;
  129.  
  130.  
  131. {       Trig Functions     }
  132.  
  133.  
  134.    FUNCTION DEGREE (X:real) : real;
  135.       BEGIN
  136.          DEGREE := X*57.295779513
  137.       END; {DEGREE}
  138.  
  139.    FUNCTION RADIAN (X:real) : real;
  140.       BEGIN
  141.          RADIAN := X/57.295779513
  142.       END; {RADIAN}
  143.  
  144.    FUNCTION SIND (X:real) : real;
  145.       BEGIN
  146.          SIND := SIN(X/57.295779513)
  147.       END; {SIND}
  148.  
  149.    FUNCTION COSD (X:real) : real;
  150.       BEGIN
  151.          COSD := COS(X/57.295779513)
  152.       END; {COSD}
  153.  
  154.    FUNCTION TAN (X:real) : real;
  155.       BEGIN
  156.          TAN := SIN(X)/COS(X)
  157.       END; {TAN}
  158.  
  159.    FUNCTION TAND (X:real) : real;
  160.       BEGIN
  161.          TAND := TAN(X/57.295779513)
  162.       END; {TAND}
  163.  
  164.    FUNCTION ARCTAND (X:real) : real;
  165.       BEGIN
  166.          ARCTAND := 57.295779513*ARCTAN(X)
  167.       END; {ARCTAND}
  168.  
  169.    FUNCTION ARCTANYX(Y,X : real) : real;
  170.       VAR
  171.          z,zmax : real;
  172.       BEGIN
  173.          z :=ArcTanD(y/x);
  174.          IF (y>0) AND (x>0) THEN zmax := 90;
  175.          IF (y>0) AND (x<0) THEN zmax := 180;
  176.          IF (y<0) AND (x<0) THEN zmax := 270;
  177.          IF (y<0) AND (x>0) THEN zmax := 360;
  178.          WHILE z>zmax DO z := z -180;
  179.          WHILE z<(zmax-90) DO z := Z + 180;
  180.          ARCTANYX := z;
  181.       END; {ARCTANYX}
  182.  
  183.    FUNCTION ARCSIN (X:real) : real;
  184.       BEGIN
  185.          IF X * X >= 1 THEN ARCSIN := 1.570796327 ELSE
  186.          ARCSIN := ARCTAN(X/SQRT(1-X*X))
  187.       END; {ARCSIN}
  188.  
  189.    FUNCTION ARCSIND (X:real) : real;
  190.       BEGIN
  191.          ARCSIND := 57.295779513 * ARCSIN(X)
  192.       END; {ARCSIND}
  193.  
  194.    FUNCTION ARCCOS (X:real) : real;
  195.       BEGIN
  196.          IF X * X >= 1 THEN ARCCOS := 1 ELSE
  197.          ARCCOS := 1.570796327-ARCTAN(X/SQRT(1-X*X))
  198.       END; {ARCCOS}
  199.  
  200.    FUNCTION ARCCOSD (X:real) : real;
  201.       BEGIN
  202.          ARCCOSD := 57.295779513*ARCCOS(X)
  203.       END;
  204.  
  205.  
  206.    PROCEDURE Frame(UpperLeftX, UpperLeftY, LowerRightX, LowerRightY: Integer);
  207.    VAR
  208.       i: Integer;
  209.    BEGIN
  210.       GotoXY(UpperLeftX, UpperLeftY);  Write('┌');
  211.       FOR i:=UpperLeftX+1 TO LowerRightX-1 DO Write('─');
  212.       Write('┐');
  213.       FOR i:=UpperLeftY+1 TO LowerRightY-1 DO
  214.       BEGIN
  215.          GotoXY(UpperLeftX , i);  Write('│');
  216.          GotoXY(LowerRightX, i);  Write('│');
  217.       END;
  218.             GotoXY(UpperLeftX, LowerRightY);
  219.       Write('└');
  220.       FOR i:=UpperLeftX+1 TO LowerRightX-1 DO Write('─');
  221.       Write('┘');
  222.    END  { Frame };
  223.  
  224.  
  225. {      Astronomy Functions     }
  226.  
  227.  
  228.    FUNCTION Minutes (time :real) : integer;
  229.        BEGIN
  230.             Minutes := ABS(TRUNC(60*(FRAC(time))))
  231.        END;
  232.  
  233.    FUNCTION Seconds (time :real) : integer;
  234.        BEGIN
  235.             Seconds := ABS(TRUNC(60*(FRAC(60*FRAC(time)))))
  236.        END;
  237.  
  238.    PROCEDURE deg_min_sec (VAR angle:real);
  239.    BEGIN
  240.       IF NOT Display_Degrees_in_Decimal
  241.          THEN
  242.             BEGIN
  243.                IF angle < 0 THEN write('-') ELSE write(' ');
  244.                write(abs(trunc(angle)):3,'° ',minutes(angle):2,''' ',
  245.                      seconds(angle):2,'"');
  246.             END
  247.          ELSE
  248.             write(angle:8:4,'° ');
  249.        END;
  250.  
  251.    PROCEDURE hours_min_sec (VAR angle:real);
  252.    BEGIN
  253.       IF NOT Display_Degrees_in_Decimal
  254.          THEN
  255.             BEGIN
  256.                IF angle < 0 THEN write('-') ELSE write(' ');
  257.                write(abs(trunc(angle/15)):3,'h ',minutes(angle/15):2,'m ',
  258.                      seconds(angle/15):2,'s');
  259.             END
  260.          ELSE
  261.             write(angle/15:8:4,'h');
  262.        END;
  263.  
  264.    FUNCTION Daylight_Savings : boolean;
  265.      BEGIN
  266.        Daylight_Savings := (((Month + Day/100) < Daylight_savings_end)  AND
  267.                            ((Month + Day/100) > Daylight_savings_start) AND
  268.                            Daylight_savings_Correction_Enabled);
  269.      END;
  270.  
  271.    FUNCTION Anomaly (M,Ecc : real) : real;
  272.        VAR
  273.           d : real;
  274.           E : real;
  275.        BEGIN
  276.           E := M;
  277.           d := E - Ecc * SIN(E) - M;
  278.           WHILE abs(d) > 0.000001 DO
  279.              BEGIN
  280.                 E := E - d/(1-Ecc * COS(E));
  281.                 d := E - Ecc * SIN(E) - M
  282.              END;
  283.           Anomaly := 2 * ArcTan(sqrt((1+Ecc)/(1-Ecc)) * Tan(E/2))
  284.         END;{Anomaly}
  285.  
  286.    FUNCTION JulianDay (month,day,year :integer) : real;
  287.        VAR
  288.           A,B,C,D : real;
  289.  
  290.        BEGIN
  291.          IF (month=1) OR (month=2)
  292.             THEN
  293.                BEGIN
  294.                   year := year-1;
  295.                   month := month+12
  296.                END;
  297.          IF (year>1582) OR (year=1582) AND (month>10) OR (year=1582) AND (month=10) AND (day>15)
  298.             THEN
  299.                BEGIN
  300.                   A := INT(year/100);
  301.                   B := 2-A+INT(A/4)
  302.                END
  303.             ELSE
  304.                BEGIN
  305.                   B := 0
  306.                END;
  307.  
  308.         C := INT(365.25*year);
  309.         D := INT(30.6001*(month+1));
  310.         JulianDay := B+C+D+day+1720994.5;
  311.       END;{JulianDay}
  312.  
  313. FUNCTION ECL_TO_RA (L,B : real) : real;
  314.    VAR
  315.      x : real;
  316.      y : real;
  317.      z : real;
  318.      zmax : real;
  319.    BEGIN
  320.      x := COSD(L);
  321.      y := SIND(L) * 0.91746406 - TAND(B) * 0.397818676;
  322.      z := ArcTanYX(y,x);
  323.      ECL_TO_RA := Z;
  324.    END;{ECL_TO_RA}
  325.  
  326. FUNCTION ECL_TO_DEC (L,B : real) : real;
  327.    BEGIN
  328.      ECL_TO_DEC := ARCSIND(SIND(B) * 0.91746406 + COSD(B) * SIND(L)
  329.                    * 0.397818676);
  330.    END;
  331.  
  332. FUNCTION LST_TO_LCT (LST,Long : real;
  333.                      Year, Day_of_year : integer;
  334.                      Zone_Corr : real) : real;
  335.    VAR
  336.      B,T,GST,GMT : real;
  337.    BEGIN
  338.         GST := LST + Long/15;
  339.         IF GST > 24 THEN GST := GST - 24;
  340.         IF GST < 0 THEN GST := GST + 24;
  341.         T := (JulianDay(1,0,year) - 2415020.0)/36525.0;
  342.         B := 24 - 6.6460656 - (2400.051262 * T) - (0.00002581 * T * T) +
  343.              (24 * (year - 1900));
  344.         T := GST -  Day_of_Year  * 0.0657098 + B;
  345.         IF T > 24 THEN T := T - 24;
  346.         IF T < 0 THEN T := T + 24;
  347.         GMT := T * 0.99727;
  348.         LST_TO_LCT := GMT - Zone_Corr;
  349.     END;
  350.  
  351. PROCEDURE make_degrees_in_range(VAR n : real);
  352. BEGIN
  353.    WHILE n > 360 DO
  354.       n := n - 360;
  355.    WHILE n < 0 DO
  356.       n := n + 360;
  357. END;
  358.  
  359. PROCEDURE Make_Hours_in_Range(VAR n : real);
  360.    BEGIN
  361.       IF n > 24 THEN n := n - 24;
  362.       IF n <  0 THEN n := n + 24;
  363.    END;
  364.  
  365. PROCEDURE time_window;
  366.    BEGIN
  367.    WINDOW(1,1,80,25);
  368.    frame(56,1,80,10);
  369.    WINDOW(58,2,78,10);
  370.    GotoXY(1,1);
  371.    writeln('    ',month,'-',day,'-',year);
  372.    JD := Julianday(month,day,year);
  373.    daylight_savings_start :=
  374.           4 + (30 - Round(7 * Frac((Julianday(4,30,year) + 1.5)/7))) * 0.01;
  375.    daylight_savings_end :=
  376.           10 + (31 - Round(7 * Frac((Julianday(10,31,year) + 1.5)/7))) * 0.01;
  377.    writeln(' JD =',JD:10:0);
  378.    write('Long = ');
  379.    deg_min_sec(Longitude);
  380.    writeln;
  381.    write('Lat  = ');
  382.    deg_min_sec(Latitude);
  383.    writeln;
  384.    LCT := Hour + Minute/60 + 1/120;
  385.    writeln('    LCT = ',Hour:2,'h ',Minute:2,'m');
  386.    GMT := Zone_correction + LCT;
  387.    IF  Daylight_Savings THEN
  388.       GMT := GMT - 1;
  389.    Make_Hours_in_Range(GMT);
  390.    writeln ('    GMT = ',trunc(GMT):2,'h ',minutes(GMT):2,'m');
  391.    T := (JulianDay(1,0,year) - 2415020.0)/36525.0;
  392.    B := 24 - 6.6460656 - (2400.051262 * T) - (0.00002581 * T * T) +
  393.         (24 * (year - 1900));
  394.    GST := 0.0657098 * Day_of_year - B + GMT * 1.002738;
  395.    Make_Hours_in_Range(GST);
  396.    writeln ('    GST = ',trunc(GST):2,'h ',minutes(GST):2,'m');
  397.    LST := GST - Longitude/15;
  398.    Make_Hours_in_Range(LST);
  399.    writeln ('    LST = ',trunc(LST):2,'h ',minutes(LST):2,'m');
  400.    END;{time_window}
  401.  
  402. PROCEDURE Locate_Position_of_Planet_in_Its_Own_Orbital_Plane(n:integer);
  403.    BEGIN
  404.     WITH planets[n] DO
  405.      BEGIN
  406.        planet_long_of_ascend_node := (360/365.2422) * (Julianday(month,day,year)
  407.                                   - Julianday(1,0,Epoch)) / period;
  408.        Make_Degrees_in_Range(planet_long_of_ascend_node);
  409.        planet_mean_anomaly := planet_long_of_ascend_node + long_at_epoch
  410.                               - long_of_per;
  411.        planet_true_anomaly := DEGREE(ANOMALY(RADIAN(planet_mean_anomaly),Ecc));
  412.        planet_heliocentric_long := planet_true_anomaly + long_of_per;
  413.        Make_Degrees_in_Range(planet_heliocentric_long);
  414.        planet_radius_vector := axis*(1-Ecc*Ecc)/(1+Ecc*COSD(planet_true_anomaly));
  415.        planet_helio_ecliptic_lat := ARCSIND(SIND(planet_heliocentric_long -
  416.                                     long_of_ascend_node) * SIND(Inc));
  417.     END;
  418.  
  419.    END;
  420.  
  421.  
  422. PROCEDURE Plot_planet(n,scale : integer);
  423. FUNCTION compute_radius_vector(n,angle : integer): real ;
  424. BEGIN
  425.    compute_radius_vector := planets[n].axis *
  426.                             (1 - planets[n].ecc * planets[n].ecc)
  427.                             /(1 + planets[n].ecc * COSD(angle
  428.                             - planets[n].long_of_per));
  429. END;
  430. VAR
  431. delX,delY,Xold,Xnew,Yold,Ynew,A : integer;
  432. BEGIN
  433.    Locate_Position_of_Planet_in_Its_Own_Orbital_Plane(n);
  434.    Ynew := Round(scale * planet_radius_vector * SIND(planet_heliocentric_long))
  435.            + CenterY;
  436.    Xnew := Round(scale * Screen_Aspect * planet_radius_vector *
  437.            COSD(planet_heliocentric_long)) + CenterX;
  438.    FOR delX := -3 TO 3 DO
  439.       FOR delY := -1 TO 1 DO
  440.           plot(Xnew + delX,Ynew + delY,1);
  441.    GoToXY((Xnew DIV 8)+2,(Ynew DIV 8)+1);
  442.    write(planets[n].name);
  443.  
  444.    Xnew := CenterX + Round(compute_radius_vector(n,0) * scale * Screen_Aspect);
  445.    Ynew := CenterY;
  446.    FOR A := 1 TO 30 DO
  447.       BEGIN
  448.          Xold := Xnew;
  449.          Yold := Ynew;
  450.          Xnew := Round(compute_radius_vector(n,A*12) * scale * Screen_Aspect
  451.                  * COSD(A * 12)) + CenterX;
  452.          Ynew := Round(compute_radius_vector(n,A*12) * scale * SIND(A * 12))
  453.                  + CenterY;
  454.          Draw(Xold,Yold,Xnew,Ynew,1);
  455.       END;
  456. END;
  457.  
  458. {    Orbital Data for Planets   }
  459.  
  460.  
  461. PROCEDURE Orbital_Data;
  462. BEGIN
  463. WITH planets[1] DO
  464.    BEGIN
  465.       name :='Mercury   ';
  466.       Epoch := 1980;
  467.       Period := 0.24085;
  468.       Long_at_Epoch := 231.2973;
  469.       Long_of_Per := 77.1442128;
  470.       Ecc := 0.2056306;
  471.       Axis := 0.3870986;
  472.       Inc := 7.0043579;
  473.       Long_of_Ascend_Node := 48.0941733;
  474.       Size := 6.74;
  475.       Brightness := 0.000001918;
  476.    END;
  477. WITH planets[2] DO
  478.    BEGIN
  479.       name :='Venus     ';
  480.       Epoch := 1980;
  481.       Period := 0.61521;
  482.       Long_at_Epoch := 355.73352;
  483.       Long_of_Per := 131.2895792;
  484.       Ecc := 0.0067826;
  485.       Axis := 0.7233316;
  486.       Inc := 3.3944359;
  487.       Long_of_Ascend_Node := 76.4997524;
  488.       Size := 16.92;
  489.       Brightness := 0.00001721;
  490.    END;
  491. WITH planets[3] DO
  492.    BEGIN
  493.       name :='Earth     ';
  494.       Epoch := 1980;
  495.       Period := 1.00004;
  496.       Long_at_Epoch := 98.833540;
  497.       Long_of_Per := 102.596403;
  498.       Ecc := 0.016718;
  499.       Axis := 1.0;
  500.       Inc := 0;
  501.       Long_of_Ascend_Node := 0;
  502.       Size := 17;
  503.       Brightness := 0;
  504.    END;
  505. WITH planets[4] DO
  506.    BEGIN
  507.       name :='Mars      ';
  508.       Epoch := 1980;
  509.       Period := 1.88089;
  510.       Long_at_Epoch := 126.30783;
  511.       Long_of_Per := 335.6908166;
  512.       Ecc := 0.0933865;
  513.       Axis := 1.5236883;
  514.       Inc := 1.8498011;
  515.       Long_of_Ascend_Node := 49.4032001;
  516.       Size := 9.36;
  517.       Brightness := 0.00000454;
  518.    END;
  519. WITH planets[5] DO
  520.    BEGIN
  521.       name :='Jupiter   ';
  522.       Epoch := 1980;
  523.       Period := 11.86224;
  524.       Long_at_Epoch := 146.966365;
  525.       Long_of_Per := 14.0095493;
  526.       Ecc := 0.0484658;
  527.       Axis := 5.202561;
  528.       Inc := 1.3041819;
  529.       Long_of_Ascend_Node := 100.2520175;
  530.       Size := 196.74;
  531.       Brightness := 0.0001994;
  532.    END;
  533. WITH planets[6] DO
  534.    BEGIN
  535.       name :='Saturn    ';
  536.       Epoch := 1980;
  537.       Period := 29.4571;
  538.       Long_at_Epoch := 165.322242;
  539.       Long_of_Per := 92.6653974;
  540.       Ecc := 0.0556155;
  541.       Axis := 9.554747;
  542.       Inc := 2.4893741;
  543.       Long_of_Ascend_Node := 113.4888341;
  544.       Size := 165.60;
  545.       Brightness := 0.0001740;
  546.    END;
  547.  
  548. WITH planets[7] DO
  549.    BEGIN
  550.       name :='Uranus    ';
  551.       Epoch := 1980;
  552.       Period := 84.01247;
  553.       Long_at_Epoch := 228.0708551;
  554.       Long_of_Per := 172.7363288;
  555.       Ecc := 0.0463232;
  556.       Axis := 19.21814;
  557.       Inc := 0.7729895;
  558.       Long_of_Ascend_Node := 73.8768642;
  559.       Size := 65.80;
  560.       Brightness := 0.00007768;
  561.    END;
  562. WITH planets[8] DO
  563.    BEGIN
  564.       name :='Neptune   ';
  565.       Epoch := 1980;
  566.       Period := 164.79558;
  567.       Long_at_Epoch := 260.3578998;
  568.       Long_of_Per := 47.8672148;
  569.       Ecc := 0.0090021;
  570.       Axis := 30.10957;
  571.       Inc := 1.7716017;
  572.       Long_of_Ascend_Node := 131.5606494;
  573.       Size := 62.20;
  574.       Brightness := 0.00007597;
  575.    END;
  576.  
  577. WITH planets[9] DO
  578.    BEGIN
  579.       name :='Pluto     ';
  580.       Epoch := 1980;
  581.       Period := 250.9;
  582.       Long_at_Epoch := 209.439;
  583.       Long_of_Per := 222.972;
  584.       Ecc := 0.25387;
  585.       Axis := 39.78459;
  586.       Inc := 17.137;
  587.       Long_of_Ascend_Node := 109.941;
  588.       Size := 8.2;
  589.       Brightness := 0.000004073;
  590.    END;
  591. END; {Orbital_Data}
  592.  
  593. PROCEDURE Show_Menu;
  594. BEGIN
  595.    WINDOW(1,1,80,25);
  596.    CLRSCR;
  597.    Time_Window;
  598.    WINDOW(1,1,80,25);
  599.    GotoXY(1,1);
  600.  
  601.    FOR num := 1 TO Number_of_Planets DO
  602.       BEGIN
  603.          WITH planets[num] DO
  604.             writeln(num,'  ',name);
  605.       END;
  606.    writeln;
  607.    writeln('D  Change Date/Time');
  608.    writeln('I  Plot Inner Planets');
  609.    writeln('L  Change Long/Lat');
  610.    writeln('M  Menu');
  611.    writeln('O  Plot Outer Planets');
  612.    writeln('S  Set_up');
  613.    writeln('Q  Quit');
  614.    writeln;
  615.    writeln('Enter command:');
  616. END; {Show_menu}
  617.  
  618. PROCEDURE Set_Up;
  619. BEGIN
  620.    writeln('Display Degrees in Decimal Degrees ( Y/N )? ');
  621.    Read(KBD,Ch);
  622.    IF (Ch = 'Y') OR (Ch = 'y')
  623.       THEN Display_Degrees_in_Decimal := true;
  624.    IF (Ch = 'N') OR (Ch = 'n')
  625.       THEN Display_Degrees_in_Decimal := false;
  626.    writeln('Correct for daylight savings ( Y/N )?');
  627.    Read(KBD,Ch);
  628.    IF (Ch = 'Y') OR (Ch = 'y')
  629.       THEN Daylight_Savings_Correction_Enabled := true;
  630.    IF (Ch = 'N') OR (Ch = 'n')
  631.       THEN Daylight_Savings_Correction_Enabled := false;
  632.    writeln('Display Hour Angle in Degrees ( Y/N )?');
  633.    Read(KBD,Ch);
  634.    IF (Ch = 'Y') OR (Ch = 'y')
  635.       THEN Display_Hour_Angle_in_Degrees := true;
  636.    IF (Ch = 'N') OR (Ch = 'n')
  637.       THEN Display_Hour_Angle_in_Degrees := false;
  638.    writeln('Display RA in Degrees ( Y/N )?');
  639.    Read(KBD,Ch);
  640.    IF (Ch = 'Y') OR (Ch = 'y')
  641.       THEN Display_RA_in_Degrees := true;
  642.    IF (Ch = 'N') OR (Ch = 'n')
  643.       THEN Display_RA_in_Degrees := false;
  644. END;
  645.  
  646.  
  647. PROCEDURE Plot_Inner_Planets;
  648. BEGIN
  649.    HiRes;
  650.    HiResColor(1);
  651.    Window(1,1,80,25);
  652.    GotoXY(62,1);
  653.    writeln('    ',month,'-',day,'-',year);
  654.    GotoXY(62,2);
  655.    writeln(' JD =',JD:10:0);
  656.    plot(CenterX,CenterY,1);
  657.    FOR num := 1 TO 4 DO
  658.       plot_planet(num,scale_for_inner_planets);
  659. END;
  660.  
  661.  
  662. PROCEDURE Plot_Outer_Planets;
  663. BEGIN
  664.    HiRes;
  665.    HiResColor(1);
  666.    Window(1,1,80,25);
  667.    GotoXY(62,1);
  668.    writeln('    ',month,'-',day,'-',year);
  669.    GotoXY(62,2);
  670.    writeln(' JD =',JD:10:0);
  671.    plot(CenterX,CenterY,1);
  672.    plot_planet(3,scale_for_outer_planets);
  673.    FOR num := 5 TO 9 DO
  674.       plot_planet(num,scale_for_outer_planets);
  675. END;
  676.  
  677.  
  678. PROCEDURE Change_Date_Time;
  679.      BEGIN
  680.        time_window;
  681.        WINDOW(1,1,80,25);
  682.        writeln;
  683.        writeln ('Enter month  day  year');
  684.        readln (month,day,year);
  685.        writeln ('Enter hour minute');
  686.        readln (hour,minute);
  687.        day_of_year := TRUNC(Julianday(month,day,year)-Julianday(1,0,year));
  688.        time_window;
  689.      END;
  690.  
  691.  
  692.  
  693. PROCEDURE Change_Long_Lat;
  694.     BEGIN
  695.       time_window;
  696.       WINDOW(1,1,80,25);
  697.       writeln;
  698.       writeln ('Enter Longitude  Latitude');
  699.       readln (longitude,Latitude);
  700.       Zone_correction := trunc(longitude/15) + 1;
  701.       writeln ('Enter zone correction (Default = ',zone_correction,' ):');
  702.       readln (zone_correction);
  703.       time_window;
  704.     END;
  705.  
  706.  
  707. PROCEDURE Locate_Planet;
  708. LABEL
  709. End_of_Locate_Planet;
  710. BEGIN
  711.    Time_window;
  712.    WINDOW(1,1,80,25);
  713.    frame(1,1,54,9);
  714.    WINDOW(3,2,53,10);
  715.    GotoXY(1,1);
  716.    Locate_Position_of_Planet_in_Its_Own_Orbital_Plane(num);
  717.    WITH planets[num] DO
  718.      BEGIN
  719.        writeln('                          ',name);
  720.        write('long of ascend node  = ');
  721.        deg_min_sec(planet_long_of_ascend_node);
  722.        writeln;
  723.        write('mean anomaly         = ');
  724.        deg_min_sec(planet_mean_anomaly);
  725.        writeln;
  726.        write('true anomaly         = ');
  727.        deg_min_sec(planet_true_anomaly);
  728.        writeln;
  729.        write('heliocentric long    = ');
  730.        deg_min_sec(planet_heliocentric_long);
  731.        writeln;
  732.        writeln('radius vector        = ',planet_radius_vector:8:3,' AU');
  733.        write('helio ecliptic lat.  =');
  734.        deg_min_sec(planet_helio_ecliptic_lat);
  735.        writeln;
  736.      END;
  737.  
  738. {  Locate Position of Earth in Its Own Orbital Plane   }
  739.  
  740.    IF num = 3 THEN GoTO End_of_Locate_Planet;
  741.    WINDOW(40,2,53,10);
  742.    GotoXY(1,1);
  743.    writeln('     Earth');
  744.    WITH planets[3] DO
  745.      BEGIN
  746.  
  747.        earth_long_of_ascend_node := (360/365.2422) * (Julianday(month,day,year)
  748.                                   - Julianday(1,0,Epoch))/period;
  749.        Make_Degrees_in_Range(earth_long_of_ascend_node);
  750.        deg_min_sec(earth_long_of_ascend_node);
  751.        writeln;
  752.  
  753.        earth_mean_anomaly := earth_long_of_ascend_node + long_at_epoch
  754.                               - long_of_per;
  755.        deg_min_sec(earth_mean_anomaly);
  756.        writeln;
  757.  
  758.        earth_true_anomaly := DEGREE(ANOMALY(RADIAN(earth_mean_anomaly),Ecc));
  759.        deg_min_sec(earth_true_anomaly);
  760.        writeln;
  761.  
  762.        earth_heliocentric_long := earth_true_anomaly + long_of_per;
  763.        Make_Degrees_in_Range(earth_heliocentric_long);
  764.        deg_min_sec(earth_heliocentric_long);
  765.        writeln;
  766.  
  767.        earth_radius_vector := axis*(1-Ecc*Ecc)/(1+Ecc*COSD(earth_true_anomaly));
  768.        writeln(earth_radius_vector:8:3,' AU');
  769.        zero := 0;
  770.        deg_min_sec(zero);
  771.      END;
  772.  
  773.  
  774. {  Project Position of Planet onto plane of ecliptic   }
  775.  
  776.  
  777.    WINDOW(1,1,80,25);
  778.    frame(1,10,54,13);
  779.    WINDOW( 12,11,79,23);
  780.    GotoXY(1,1);
  781.    WITH planets[num] DO
  782.      BEGIN
  783.         y := SIND(planet_heliocentric_long - long_of_ascend_node) * COSD(Inc);
  784.         x := COSD(planet_heliocentric_long - long_of_ascend_node);
  785.         z :=ArcTanYX(y,x);
  786.         projected_heliocentric_long := Z + long_of_ascend_node;
  787.         Make_Degrees_in_Range(projected_heliocentric_long);
  788.         write('Projected Longitude  = ');
  789.         deg_min_sec(projected_heliocentric_long);
  790.         writeln;
  791.         Projected_Radius_vector := planet_radius_vector *
  792.                                    COSD(planet_helio_ecliptic_lat);
  793.         writeln('Projected Radius     = ',Projected_Radius_Vector:8:3,' AU');
  794.       END;
  795.  
  796. {    Calculate Ecliptical Coordinates   }
  797.  
  798.  
  799.    WINDOW(1,1,80,25);
  800.    frame(1,18,25,24);
  801.    WINDOW(3,19,24,23);
  802.    GotoXY(1,1);
  803.    IF (Planet_Radius_Vector < Earth_Radius_Vector)
  804.    THEN
  805.       Geocentric_Ecliptic_Long := Earth_Heliocentric_long + 180 +
  806.            ArcTanD(Projected_Radius_Vector * SIND(Earth_Heliocentric_long -
  807.            Projected_heliocentric_long)/(Earth_Radius_Vector -
  808.            Projected_Radius_Vector * COSD(Earth_Heliocentric_long -
  809.            Projected_heliocentric_long)))
  810.    ELSE
  811.       Geocentric_Ecliptic_Long := Projected_Heliocentric_Long +
  812.            ArcTanD(Earth_Radius_Vector * SIND(Projected_heliocentric_long -
  813.            Earth_Heliocentric_Long) / (Projected_Radius_Vector -
  814.            Earth_Radius_Vector * COSD(Projected_heliocentric_long -
  815.            Earth_Heliocentric_Long)));
  816.     Make_Degrees_in_Range(Geocentric_Ecliptic_Long);
  817.    writeln('      ECLIPTIC');
  818.    writeln;
  819.    write(' Long = ');
  820.    deg_min_sec(Geocentric_Ecliptic_Long);
  821.    writeln;
  822.    Geocentric_latitude := ArcTanD(Projected_Radius_Vector *
  823.                           TAND(planet_helio_ecliptic_lat) *
  824.                           SIND(Geocentric_Ecliptic_Long -
  825.                           Projected_Heliocentric_Long) / (Earth_Radius_Vector *
  826.                           SIND(Projected_Heliocentric_Long -
  827.                           Earth_Heliocentric_Long)));
  828.    write(' Lat  = ');
  829.    deg_min_sec(Geocentric_latitude);
  830.    writeln;
  831.  
  832.  
  833. {     Calculate Equatorial Coordinates    }
  834.  
  835.    WINDOW(1,1,80,25);
  836.    frame(27,18,52,24);
  837.    WINDOW(29,19,51,24);
  838.    GotoXY(1,1);
  839.    writeln('     EQUATORIAL');
  840.    writeln;
  841.    RA := ECL_TO_RA(Geocentric_Ecliptic_Long,Geocentric_Latitude);
  842.    write(' RA  = ');
  843.    IF Display_RA_in_Degrees
  844.    THEN
  845.       deg_min_sec(RA)
  846.    ELSE
  847.       hours_min_sec(RA);
  848.    writeln;
  849.    DEC := ECL_TO_DEC(Geocentric_Ecliptic_Long,Geocentric_Latitude);
  850.    write(' DEC = ');
  851.    deg_min_sec(DEC);
  852.    writeln;
  853.    Hour_Angle := LST - RA/15;
  854.    IF Hour_Angle < 0 THEN Hour_Angle := Hour_Angle + 24;
  855.    Hour_Angle_in_Degrees := Hour_Angle * 15;
  856.    write(' HA  = ');
  857.    IF Display_Hour_Angle_in_Degrees
  858.    THEN
  859.       deg_min_sec(Hour_Angle_in_Degrees)
  860.    ELSE
  861.       hours_min_sec(Hour_Angle_in_Degrees);
  862.  
  863. {    Calculate Horizontal Coordinates   }
  864.  
  865.    WINDOW(1,1,80,25);
  866.    frame(54,18,80,24);
  867.    WINDOW(56,19,78,23);
  868.    GotoXY(1,1);
  869.    writeln('      HORIZONTAL');
  870.    writeln;
  871.    Altitude := ArcSinD(SIND(Dec) * SIND(Latitude) +
  872.                COSD(Dec) * COSD(Latitude) * COSD(Hour_Angle * 15));
  873.    write(' Alt  = ');
  874.    deg_min_sec(Altitude);
  875.    writeln;
  876.  
  877.    Azimuth := ArcCosD((SIND(Dec) - SIND(Latitude) * SIND(Altitude))/
  878.               (COSD(Latitude) * COSD(Altitude)));
  879.    IF SIND(Hour_Angle) > 0 THEN
  880.    Azimuth := 360 - Azimuth;
  881.    write(' Azim = ');
  882.    deg_min_sec(Azimuth);
  883.  
  884.  
  885. {   Calculate Time and Position of Rise and Set   }
  886.  
  887.  
  888.    WINDOW(1,1,80,25);
  889.    frame(1,14,54,17);
  890.    WINDOW(4,15,52,17);
  891.    GotoXY(1,1);
  892.    Azimuth_Rise := SIND(DEC)/COSD(Latitude);
  893.    IF (Azimuth_Rise < -1) OR (Azimuth_Rise > 1) THEN
  894.       IF Altitude < 0
  895.          THEN
  896.             writeln('never rises above horizon')
  897.          ELSE
  898.             writeln('never sets below horizon')
  899.    ELSE
  900.       BEGIN
  901.         Azimuth_Rise := ARCCOSD(Azimuth_Rise);
  902.         Azimuth_Set := 360 - Azimuth_Rise;
  903.         LST_Rise := 24 + RA/15 - (ARCCOSD(-TAND(Latitude) * TAND(DEC)))/15;
  904.         IF LST_Rise > 24 THEN LST_Rise := LST_Rise - 24;
  905.         LST_SET := RA/15 + (ARCCOSD(-TAND(Latitude) * TAND(DEC)))/15;
  906.         IF LST_SET > 24 THEN LST_SET := LST_SET - 24;
  907.  
  908.         LCT_SET := LST_TO_LCT(LST_SET,Longitude,Year,Day_of_Year,
  909.                     Zone_Correction);
  910.         IF Daylight_Savings THEN
  911.            LCT_SET := LCT_SET + 1;
  912.         IF LCT_SET < 0 THEN LCT_SET := LCT_SET + 24;
  913.  
  914.         LCT_RISE := LST_TO_LCT(LST_RISE,Longitude,Year,Day_of_Year,
  915.                     Zone_Correction);
  916.         IF Daylight_Savings THEN
  917.            LCT_RISE := LCT_RISE + 1;
  918.         IF LCT_RISE < 0 THEN LCT_RISE := LCT_RISE + 24;
  919.  
  920.         write('Rises at ',trunc(LCT_Rise):2,':',minutes(LCT_Rise):2,'  LCT');
  921.         write('      Azimuth  ');
  922.         deg_min_sec(Azimuth_Rise);
  923.         writeln;
  924.         write('Sets at  ',trunc(LCT_Set):2,':',minutes(LCT_Set):2,'  LCT');
  925.         write('      Azimuth  ');
  926.         deg_min_sec(Azimuth_Set);
  927.       END;
  928.  
  929.  
  930. {    Calculate Phase,Distance, Diameter, Magnitude   }
  931.  
  932.  
  933.    Window(1,1,80,25);
  934.    frame(56,11,80,17);
  935.    Window(57,12,79,16);
  936.    GotoXY(1,1);
  937.  
  938.    Phase := (1+COSD(Geocentric_ecliptic_long - planet_heliocentric_long))/2;
  939.    writeln('Phase     = ',100 * Phase:6:2,'%');
  940.    Distance_from_Earth := sqrt(Earth_Radius_Vector * Earth_Radius_Vector
  941.                           + Planet_Radius_Vector * Planet_Radius_Vector
  942.                           - 2 * Earth_Radius_Vector * Planet_Radius_Vector *
  943.                           COSD(planet_heliocentric_long -
  944.                           earth_heliocentric_long));
  945.    writeln('Distance  = ',Distance_from_Earth:6:2,' AU');
  946.  
  947.    WITH Planets[num] DO
  948.    Diameter := size / Distance_from_Earth;
  949.    writeln('Diameter  = ',Diameter:6:2,'"');
  950.  
  951.    WITH Planets[num] DO
  952.    Magnitude := 2.17147 * ln(Distance_from_Earth * Planet_Radius_Vector /
  953.                 Brightness * sqrt(Phase)) - 26.7;
  954.    writeln('Magnitude = ',Magnitude:6:2);
  955. End_of_Locate_Planet:
  956. END; {Locate_Planet}
  957.  
  958. BEGIN {Planets}
  959.  
  960. {      Initalize Default Parameters   }
  961.  
  962.    CLRSCR;
  963.    MSDOS_time.AX := $2C00;{time call to MSDOS}
  964.    MsDos(MSDOS_time);
  965.    Hour :=Hi(MSDOS_time.CX);
  966.    Minute := Lo(MSDOS_time.CX);
  967.    MSDOS_time.AX := $2A00; {Date call to MSDOS}
  968.    MsDos(MSDOS_time);
  969.    Year := MSDOS_time.CX;
  970.    Month := Hi(MSDOS_time.DX);
  971.    Day := Lo(MSDOS_time.DX);
  972.    Orbital_Data;
  973.    day_of_year := TRUNC(Julianday(month,day,year)-Julianday(1,0,year));
  974.  
  975. {     Main Program     }
  976.  
  977.  
  978. Show_Menu;
  979. REPEAT
  980.    Read(KBD, Ch);
  981.    num := ord(Ch) - ord('0');
  982.    TextMode;
  983.    Bios_Display_Mode := $2D;
  984.    ClrScr;
  985.    CASE Ch OF
  986.       'I','i' :  Plot_Inner_Planets;
  987.       'O','o' :  Plot_Outer_Planets;
  988.       'D','d' :  Change_Date_Time;
  989.       'L','l' :  Change_Long_Lat;
  990.       'M','m' :  Show_Menu;
  991.       'S','s' :  Set_Up;
  992.       '1','2','3','4','5','6','7','8','9'  : Locate_Planet;
  993.    END; {case}
  994.    Window(1,1,80,25);
  995.    GotoXY(1,25);
  996.    write('(Q)uit (M)enu (I)nner (O)uter (D)ate (L)ong   (1-9) Planet');
  997. UNTIL (Ch = 'Q') OR (Ch = 'q');
  998.    TextMode;
  999.    Bios_Display_Mode := $2D;
  1000.    ClrScr;
  1001. END.